Particle linear theory on a self-gravitating perturbed cubic Bravais lattice 
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Abstract 

Discreteness effects are a source of uncontrolled systematic errors of N-body simulations, which are 
used to compute the evolution of a self-gravitating fluid. We have already developed the so-called 
"Particle Linear Theory" (PLT), which describes the evolution of the position of self-gravitating 
particles located on a perturbed simple cubic lattice. It is the discrete analogue of the well-known 
(Lagrangian) linear theory of a self-gravitating fluid. Comparing both theories permits to quantify 
precisely discreteness effects in the linear regime. It is useful to develop the PLT also for other 
perturbed lattices because they represent different discretizations of the same continuous system. 
In this paper we detail how to implement the PLT for perturbed cubic Bravais lattices (simple, 
body and face-centered) in a cubic simulation box. As an application, we will study the discreteness 
effects — in the linear regime — of N-body simulations for which initial conditions have been set-up 
using these different lattices. 

PACS numbers: 98.80.-k, 05.70.-a, 02.50.-r, 05.40.-a 



I. INTRODUCTION 

An important problem in cosmology is the formation 
of the large scale structure. The key process involved 
is the gravitational clustering of coUisionless dark mat- 
ter, which is considered to be well described as a self- 
gravitating fluid for a wide range of scales (e.g. [1|). The 
complexity of these fluid equations (coupled with grav- 
ity) makes impossible to compute an analytical solution. 
There are therefore two common approaches to attack 
the problem: (i) a perturbative expansion in the density 
contrast 5{t) — p(r)/po— 1 (where p(r) is the local density 
and po its space average) , valid only at early times (or for 
scales in which the density contrast averaged over such 
scales is smaller than one) and (ii) N-body simulation, 
in which the fluid is discretized into particles (N-bodies) 
and then the evolution of the system computed applying 
simple gravity. 

N-body simulations are used to compute the evolution 
in the highly non-linear regime. A basic problem of this 
method is that there is no theory on the discreteness ef- 
fects due to the use of a finite number N of particles 
(e.g. BSIISIEH)- Generally, tests varying iV shows 
a "convergence" of the simulations. However, it is diffi- 
cult to infer how well this convergence has been achieved 
because of the lack of framework to refer to. For exam- 
ple, it is not known the dependence of the discreteness 
error with N . If the convergence is slow, numerical tests 
can indeed appear to converge when actually convergence 
has not been achieved (see e.g. 0,3). 

In Q and subsequently [10, [ill, [l3| we have started 
to develop a program to precisely fill this gap. We have 
developed a framework which allows us to calculate the 
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evolution — in the linear regime — of a systenr of self- 
interacting particles. This is the discrete counterpart of 
the well-known Fluid Linear Theory (hereafter FLT), and 
we called it Particle Linear Theory (hereafter PLT). We 
have shown that the fluid limit of the PLT is well de- 
fined and indeed it is the FLT. We have shown also how 
to quantify, in an essentially analytic way, the discrete- 
ness effects, with arbitrarily large precision. Moreover, 
availability of analytical results permits to evaluate the 
discreteness effects in the limit of infinite realizations. It 
avoids, in the computation of statistical quantities, the 
use of any statistical estimator and thus its subsequent 
and problematic noise. One of our conclusions was that, 
for the set-up of the initial conditions (IC), the body 
centered cubic lattice could be a better choice than the 
simple cubic (sc) one, because it might produce less dis- 
creteness effects. We will see in this paper that it is 
indeed the case in this context of linear theory. 

Moreover, another important motivation of this pa- 
per is the study of the discreteness effects in the nan — 
perturbative regime. In the forthcoming paper [13J we 
sample the same continuous field using different lattices 
and then evolve them using N-body simulations. The 
differences between the result of these simulations give 
an estimate of the lower bound of the discreteness effects 
in the non-perturhative regime. Because these differences 
are small — typically of order of a few percent in the 
power spectrum for times and scales relevant to cosmo- 
logical simulations — , an implementation of the PLT for 
these lattices is an essential tool to check that these differ- 
ences are actually discreteness effects and not numerical 
errors, finite-size effects, estimator-related errors, etc. 

In this paper we present the PLT method applied to 
any cubic Bravais lattice, i.e., to a simple cubic (sc), 
body-centered (bcc) and face-centered (fee) lattices. In 
the first section we give an summary of the PLT. Further 
details can be found in [fOll . In the following section, we 



explicitly give the details of the PLT for a sc, bee and foe 
lattices. We use Fast Fourier Transform techniques, in a 
cubic box, which is a-priori non trivial. In the last section, 
we present some applications of the method, comparing 
discreteness effects using a perturbed sc, bcc or fee lat- 
tice to set-up the IC. It is a generalization to a bcc and 
fee lattices of the work presented in [IJl . 



II. LINEARIZATION OF GRAVITY ON A 
PERTURBED LATTICE 

In this section we present a summary of the general 
method we have developed in [10, \1M to calculate the 
evolution of self-gravitating particles perturbed off a per- 
fect lattice. 

Let us consider a parallelepiped of volume V with N 
lattice sites, which are generated combining linearly the 
three primitive lattice vectors ai, a2 and a.^: 

R = R(ni,n2,n3) = ^(niai + n2a.2 -l-nsaa), (1) 



where 



G [0, Ni -i]nz 



(2) 



and i is the typical "lattice spacing"^ (we have chosen 
SLi to be dimensionless) . The total number of particles in 
the system is TV = N1N2N3 and the box a parallelepiped 
with sides A = {Ai,A2, A3}, where A^ = TVia^. 

We perform a displacement of the particles about their 
lattice position R and we write their new position r{t) 
as: 



r(t) =R + u(R,t), 



(3) 



which will evolve under the effect of gravity and where 
u(R, t) is a displacement field evaluated at the lattice 
positions. 



Definition and linearization of the gravitational 
force 



In order to have a translationally invariant system^ we 
take periodic boundary conditions. We use the method 
of replicas to compute the gravitational force. It consists 
in calculating the force not only considering the particles 
in the box of volume V but also all its images, i.e., an in- 
finite number of copies of the system. This is a standard 
scheme in cosmological N-body simulations to evaluate 
the force (see e.g. [lj|). For a well defined gravitational 



^ For a sc lattice £ is the actual lattice spacing while this is not 
true in the bcc and fee case, because all the lattice sites are not 
at the same distance each one from another. 

^ This is not to have any privileged point in the system. 



force in the infinite volume limit, it is necessary to in- 
troduce a neutralizing background which, in cosmology, 
is naturally introduced in the context of an expanding 
universe (see e.g. [V\). 

The gravitational force is linearized by expanding in 
Taylor series at linear order in the variable u(R, t) about 
the lattice position R (for more details see e.g. 1(j1). 
It is convenient to use of the dynamical matrix 2?(R) to 
express the linearized force: 



F(r) = Vl?(R-R')u(R'). 



R 



(4) 



The expression of the dynamical matrix for a generic in- 
teraction potential v{r) is USJ^: 



P^,(R^0) = c»^9,u.(R) 



where 



P^^(R = 0) 



df,d^w{ro) 



R'#0 



d^ w(v) 



9^a,ii;(R') 



dr^,dr^ 



and wij) is the periodic function defined as 
wiv) — \J t;(r + n • A), 



(5a) 
(5b) 



(6) 



(7) 



i.e., the potential due to a single particle and all its copies. 
For the gravity force, we have u(r) = —Gm/r and Eq. ([7]) 
is implicitly understood to be regularized by the addition 
of a uniform negative background. However, the sum 
(l7|) is numerically slowly convergent (it is necessary to 
sum over a huge number of replicas). To speed-up the 
computation we use the standard Ewald method, which 
consists in dividing the sum in a short range part and a 
long range one introducing a damping function C: 



(r) = V" t;(r + n • A)C(|r + n • A|, a) 



^«(r + n.A)[l-C(|r + n.A|,a)], 



(8) 



where a is a damping parameter from which the result is 
independent. A common choice for a 1/r potential is 

C(|r|,a) =erfc(a|r|). (9) 

The expression for the function w is then: 

w(r)=wM(r)-FwW(r) (10) 

and 

'('■^r) ^-GmS^ i -erfc(a|r + n-A|), (11a) 

n ' 

47r ■, 



W"' 



.W(r) = -Gm^5:^exp('-g')cos[k.r] 



4q2 



(lib) 



The Fourier vectors k are generated combining linearly 
the (dimensionless) primitive vectors in reciprocal space 
b, 



k£ = mi h m2 

where rrii are integers and 



^2 



ms 



N. 



Sii ■ hj = 2TTSij 



(12) 



(13) 



{6ij is the Kronecker delta). 
quency as 



We define the Nyquist fre- 






(14) 



It is simple to show (e.g [ifll) that the term k = is 
not included in the sum (jllbp due to the presence of 
the neutralizing background (or the space expansion in 
the cosmological context). An explicit expression of the 
dynamical matrix is given in App. VK\ 



We denote u(k, t) as the Fourier transform (hereafter FT) 
on the lattice of u(R, t) 



u(k,i)=^u(R,t)e-^''-^, 



(17) 



R 



where the sum is restricted to the simulation box (i.e. 
without considering the replicas). Using Eqs. (fT5| and 
(fT6)l we obtain the 3x3 eigenvalue problem 



fi(k,i) =2?(k)u(k,t), 



(18) 



where 2?(k) is defined analogously to u(k, t). We can 
easily diagonalize (numerically) Eq. (fTS]) . obtaining for 
each k the eigenvalue equation 



P(k)e„(k)=47rGpo£(k)e„(k), 



(19) 



where po is the average mass density po — n/V and 
£(k) the normalized eigenvalues of the dynamical matrix 
2?(k). We can decompose each mode u(k, t) in the basis 
{e„(k), n = 1,2,3} as 



B. Dynamical equations 

For simplicity we will consider a matter-dominated 
universe with zero cosmological constant (Einstein- 
deSitter, hereafter EdS)"^. This is a very good approxima- 
tion for the currently most favored ACDM cosmological 
model for the times in which PLT is a good approxima- 
tion (i.e. before shell-crossing), considering the typical 
red-shifts in which the simulations are started. The evo- 
lution of the displacement field u(R, t) is given by the 
equation 



u(k,t) = ^e„(k)/„(k,i). 



(20) 



Using Eqs. dUl), HI]) and (gOl) we get the following equa- 
tion for the coefficients /„(k, t): 

/;(k,t) + 2^/„(k,t) = l^%£M/„(k,,). (21) 

Depending on the sign of £(k), we obtain two classes of 
solutions [/„(k, t) and 14,(k, t), which are given in App. 



N 

u(R, t) = -2-u(R, t) -F ^ V P(R- R')u(R', t), (15) 



R' 



where a{t) is the scale factor and the (double) dots mean 
(double) derivative with respect to time. From Bloch 
theorem it is possible to diagonalize Eq. (fTS]) in real space 
using the following combination of plane waves: 



u(R,i) = l^u(k,i)e^'^-^ 



(16) 



where the sum is restricted to the first Brillouin zone 
(hereafter FBZ), i.e., by the set of the N vectors k'' with 
smaller modulus. These symmetrically lie around k = 0'^. 



^ For a static non-expanding universe see [Toll . 

* It is simple to show (e.g. llSlI ) that a periodic lattice with A'^ 

particles has A'^ associated independent vectors k. 
^ The FBZ is not in general symmetric about k = but this is the 

case for a cubic Bravais lattice because of the symmetries of the 

lattice. 



C. Evolution of the po^ver spectrum 

Usually, we are not interested in the position of each 
particle but in some global statistical quantities. In this 
paper, we will focus on the power spectrum (hereafter 
PS), defined as 



F(k) = hm 



V 



(22) 



where (5/5(k) is the FT of the density contrast 5p{r) = 
p{r) — po (we assume statistical homogeneity). It is pos- 
sible to show that for a small value of the displacement 
|u(R, t)\<^ i, the PS of a perturbed lattice can be writ- 
ten as ji, [li 



P(k, i) « kf,kt,gf,^{k,t), 



where 



9f^i^(y) = ,lim 



(u^(k)u:(k)) 

V 



(23) 



(24) 



Setting-up the IC at i = to in the canonical way using the 
Zcldovich approximation is equivalent to set (e.g. [9|) 

g^^(k, to) = k^ky.g(k, to). (25) 

Using Eqs. dHSJ, jHH, ^ and ^ we get: 

F(k,t)«A|,(k,t)F(k,to), (26) 

where 

Ap{k, t)=Y, K'^uA^AK t) (27) 

and (for an EdS universe) [13] 



^,,^(k,t) = ^ 



C/„(k,t) + — K(k,t) 
dro 



(,*^n j/i V*^n ji^- 



(28) 



III. DIAGONALIZATION OF THE 
DYNAMICAL MATRIX 

In this section we describe step-by-step how to diago- 
nalize the dynamical matrix. 





ai 


32 


as 


sc 


[1,0,0] 


[0,1,0]] 


[0,0,1] 


bcc 


[1,0,0] 


[0,1,0]] 


1 [1,1,1] 


fee 


i [0,1,1] 


1 [1,0,1]] 


[0,1,0] 



TABLE I: Lattiee vectors for the different kind of eubic Bra- 
vais lattiees. 





Ni 


N2 


Nb 


se 


^1/3 


^1/3 


N^/B 


bee 


(• JV\l/3 


(.)V3 


2(^)1/3 


fee 


i^r 


i^r' 


4(^)1/3 



TABLE IL Associated number particles with the lattiees vec- 
tors listed in Table |T] for the different kind of eubic Bravais 
lattiees. 



in each direction (compatible with the total number of 
particles) . We give in Table U a set of primitive lattice 
vectors and in Table |TT] the particle number associated 
with them (for a total of N particles) for a sc, bcc and 
fee lattices which fulfill the above requirements. 



A. Generation of the real space lattice 

In general, N-body simulations are performed in a cu- 
bic box, using a perturbed lattice as initial conditions. 
Therefore, to fill the simulation box in an uniform way, 
the number of particles cannot be arbitrary. In the case 
of a sc lattice, the number of points should be iV = N^^ 
(with Nsc an integer), for a bcc one N = (A'^bcc/2)^ and 
for an fee one N = (iVfcc/4)'^ (where iVbcc and iVfcc are 
also integers). 

Note that the real space vectors R, generated using 
Eq. IJ]), lie, in general, in a parallelepiped box, with sides 
{Ai, A2, A3}. Note that it is necessary to generate the 
real space vectors in this way [i.e. using Eq. ^ and ^] 
in order to use the technique of Fast Fourier Transform 
(FFT) as we will see in section UlI CI We have therefore 
to translate the R vectors into a cube using a operation 
which leaves unchanged the dynamics of the system. It 
is simple to show that the equation of motion (|15[) is 
invariant under the transformation 



R 



R 



3 



Generation of vectors in reciprocal space in the 
FBZ 



Given the primitive lattice vectors a^, the primitive 
vectors in reciprocal space are univocally defined by 
Eq. P^ . The basis we have used to generate the lat- 
tices is given in Table H] and the corresponding primitive 
reciprocal vectors are listed in Table IIIII The reciprocal 
vectors are generated using Eq. p^ where rm are the 
same integers as the ones used to generate the R vectors, 
i.e., 



G [0, Ni -i]n z. 



(30) 



^'A, 



(29) 



It is necessary, in order to use FFT techniques, to gen- 
erate the reciprocal vectors in this way, as we will see in 
section IIII CI 

However, all the k vectors used in the computation 
of the evolution of the particle position must lie in the 
FBZ (see section IITB|) but, in general, those generated 
using Eqs. p^ and (|30p do not. We can translate the 
reciprocal vectors into the FBZ using the transformation 



(where n'^ are integers). We can, then, choose three prim- 
itive lattice vectors {ai, i — 1,2,3} and the number of 
particles Ni associated with each primitive lattice vector 
(compatible with the total number of particles) which, 
using Eq. (|^^ . translate all the lattice sites into a cube. 
This is not trivial and does not work for any combina- 
tion of primitive lattice vectors and number of particles 





bi 


ba 


b3 


se 


27r[l,0,0] 


27r[0,l,0]] 


27r[0,0,l] 


bcc 


2^[1,0,-1] 


2^[0,1,-1]] 


4^[0,0,1] 


fee 


47r[-l,0, 1] 


47r[l,0,0]] 


27r[l,l,-l] 



TABLE IIL Reciprocal vectors for the different lattiees. 



sc 


bcc 


fee 


f [±1,0,0] 


f[0,+l,±l] 


^[±1,0,0] 


f[0,±l,0] 


f[+l,0,±l] 


^[0,±1,0] 


f[0,0,±l] 


f [±1,+1,0] 


2r[0,0,±l] 
f[+l, +!,+!] 


2 + 2 + 2 vectors 


4 + 4 + 4 vectors 


2 + 2 + 2 + 8 vectors 



Using Eqs. ([T]), ^^ and ^3^ we can write Eq. ^^ as 



TABLE IV: Normal vectors which define the FBZ of the bcc 
and fee lattices. 



which leads Eq. (|T8|) invariant 



k + ^ m'^i 



(31) 



where m[ are some appropriate integers. 

One can obtain a complete set of A'' k vectors which are 
in the FBZ, in the following way: compute a set of can- 
didate vectors to lie in the FBZ with Eqs. (HH) and (pij) . 
To select those which are in the FBZ, it is not efficient 
to consider the N vectors with smaller modulus because 
it is an 0{N'^) operation. The computation time for this 
can be prohibitive for large N. It is much better to con- 
struct geometrically the shape of the FBZ by considering 
some point of the reciprocal space (namely b = 0) and 
then drawing the perpendicular bisector planes of the 
translation vectors from the chosen center to the nearest 
sites of the reciprocal lattice. In Table IIVI we give the 
normal vector of this plane, with modulus equal to their 
closest distance to the center k = 0. The FBZ of the sc 
lattice is a cube of side 2n/i, the one of the bcc lattice 
a rhombic dodecahedron and the one of the fee lattice 
a cuboctahedron. Then, we select the k vectors which 
are enclosed between these planes. This is an essentially 
0{N) operation. 



C. Fast Fourier Transform 

In this section, we will carry out the EFT of some 
quantity defined on the lattice as, e.g., the dynamical 
matrix 



P(k)=^P(R) 



ik-R 



(32) 



P. 



E^" 



exp 



27ri 



niTTli n2TO2 



Ni 



N, 



^3 / 



R 



(33) 
where the indices n and m labels the R and k vectors re- 
spectively. These are the same triplets of integers which 
have been used in Eqs. dJ) and (fT^ respectively. Note 
that Eq. p3|) is a three-dimensional FT, i.e., three em- 
bedded one-dimensional FT as the one of Eq. (jCip , with 
the same running of indices [see Eqs. ^ and ([50]) ]. It is 
then straightforward to compute the FT ([55]) using any 
standard EFT routine. Note that each R vector should 
be associated in Eq. ([33]) with the indices [ni, n2, n^] with 
which it has been generated using Eq. ([1} , and not those 
that would correspond to their actual position in the cu- 
bic box after being applied the transformation ((29|l . The 
same observation holds for the k vectors, whose indices 
correspond to those used generating them with Eq. p^ . 
There exists a great number of publicly available very 
competitive EFT routines. We have used the Fastest 
Fourier Transform in the West (FFTW) [17| . since it 
can be used for any number of particles (and not only 
powers of two). 



D. Spectrum of eigenvalues of an sc, bcc and fee 
lattices 



As an application of these techniques, we show in Fig.[T] 
the spectrum of eigenvalues corresponding to a sc, bcc 
and fee lattice. We use a cubic box and Ngc — 64 for 
the sc lattice. For a comparison with the other lattices, 
we need A^bcc = 51 and iVfcc — 40. The three lattice 
presents the same branch structure (for further discus- 
sion see [101^: (i) an optical branch, with eigenvalues 
e(k^ ^- 0) = 1 and eigenvector polarized parallel to k 
(in the same limit) and (ii) two acoustic branches with 
normalized eigenvalues e(]i£ — *■ 0) = and polarized in 
the plane transverse to k (in the same limit). We also 
see that, as anticipated in [lO], the spectrum of the bcc 
and fee lattices does not present negative nor eigenvalues 
with £(k) > 1. 



IV. DISCRETENESS EFFECTS IN A BCC, FCC 
AND SC LATTICE 



where R is restricted to the simulation box. Equation 
(l32)l involves an 0{N'^) operations (an iV-term sum for 
each of the TV k vectors). However, using the so-called 
Fast Fourier Transform (FFT) technique, it is possible to 
reduce the number of operations — exploiting the sym- 
metries of the FT — to only N ln2 N operations. We give 
a brief summary of how the FFT works in App. [C) By 
using it, we can speed-up greatly the computation of the 
FTs of the dynamical matrix and the displacement field. 



In this section, we will apply the method described 
above to compare the discreteness effects when using dif- 
ferent lattices to set-up the IC, i.e., different discretiza- 
tions of the same initial density field. To do that, we 
compare FLT with PLT for the three lattices considered. 
We will study two different effects: the change in the 
amplification of the PS and the breaking of isotropy. A 
more detailed study of discreteness effects in the linear 
regime for a sc lattice can be found in [12i |. 
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FIG. 1: Normalized spectrum of eigenvalues of (from top to 
bottom) a sc, bcc and fee in function of the waveetor nor- 
malized to the Nyquist frequency of the sc lattice defined in 
Eq. (|14|) . We have performed a sampling taking 1% of the 
points. 



In the fluid limit the evolution of the PS is given by 
the well-known FLT (e.g. [H [13]): 



P""'^(fc,t) = a^(t)P(fc,io) 



(34) 




0.25 0.5 0.75 1 1.25 1.5 1.75 
k/kN 



FIG. 2: Amplification of the PS averaged over bins normalized 
to the fluid (FLT) amplification. The wavevectors have been 
normalized to the Nyquist frequency of the sc lattice (see the 
text for more details). 



evolution using PLT, is given by Eqs. ^ and (QT]). We 
set up the IC by using the ZA approximation. 

To characterize the effect of the discreteness in the am- 
plification of the PS we define the quantity 



PspiKt) 






(35) 



which is the amplification of the PS calculated with PLT 
normalized by FLT. In Fig. [2] we show the amplification 
of the PS at a{t) = 5 predicted by PLT, normahzed to the 
fluid amplification. We have averaged over 60 bins cen- 
tered in |k| with amplitude |k| ± |Ak| with |Ak| w 0.92. 
We see that the evolution of the sc lattice is slightly closer 
to the fluid one (Ap(k, i)/a(i)^ = I) for k < k^ than the 
bcc and fee lattices. This is not surprising, looking at 
the form of the spectrum of eigenvalues of the different 
lattices shown in Fig. [TJ The amplification of each sin- 
gle mode for fc <C fc^r of the PS is related essentially 
with the shape of the optical branch of the spectrum of 
eigenvalues. There are some eigenvalues in the sc lattice 
with e(k) > I which compensate, averaging over bins of 
same |k|, the largest part of eigenvalues with e(k) < 1. 
However, for k > kpf, the evolution of the sc lattice is 
farther than the other two from the fluid evolution. In 
fact, the modes with e(k) > 1 do not exist anymore for 
k > k^. Therefore, looking at the amplification of the 
PS, we can say that the sc lattice is slightly closer to the 
fluid evolution for k < fc^r. However, as we will see be- 
low, the anisotropy introduced by the sc lattice is much 
larger than the one introduced by the bcc or fee lattice. 
Let us consider the normalized dispersion of the am- 
plification of the PS, defined as in [I4I 



[we consider an initial PS which is statistically isotropic 
and we have used that a{to) = 1, see Eq. (|B2I) ]. The 



APs,{Kt)=(^M^, 



2 \ 1/2 

Psp (fc,t)\ 



Psp {k,t) 



(36) 
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FIG. 3: Normalized dispersion of the PS amplification 
^Psp{k,t) averaged over bins of the same modulus of k. 



FIG. 4: Deviation of the polarization of the eigenvectors from 
the fluid limit {6aniso(k,t > io))- 



where the average, of any function X(k, t) on the recip- 
rocal lattice, is defined as 



nk,t) = ^ J2 x{k,t), 



Nk 



(37) 



k.|k|=fc 



Nk being the number of eigenmodes at a given k. This 
quantity gives a measure of the anisotropy of the PS am- 
plification. In a system which respects the isotropy of 
a fluid, the amplification of plane waves with the same 
corresponding wavevector, but different direction, should 
be the same. In Fig. [3] we show that the sc presents a 
dispersion which is about an order of magnitude larger 
than the one of the bcc or fee, which is very similar. The 
behavior of the dispersion APsp{k,t) ex k'^ as predicted 
by PLT because the eigenvalues, for |k| < k^, goes as 
e(k) ~ 1 — a{ii)k'^/kjf, where a(k) is a function which 
depends on the particular lattice (for more details see 

Another method to quantify explicitly the breaking of 
isotropy consists in measuring the deviation of the eigen- 
values corresponding to the optical branch with the po- 
larization of them in the fiuid limit in the direction k (see 
|10lll3|). We can quantify it using the expression: 



^aniso(.k, t) — 



|u(k,t)-k-u(k,t)k|' 
|k-ii(k,i)P 



(38) 



In the infinite realizations limit, assuming that the IC 
have been set-up using the ZA [i.e. the expression 
holds], using Eq. (p8| we have: 



(D,, 



.(k,i)) 



ZA 



A A k \<: 



(-4^1/k^ky j 



(39) 



For sufficiently long times (i.e 

"^dyn = l/V47rG'po) the expression 

on how the IC have been set-up and on the cosmologi- 

cal model. It depends only on the eigenvectors, i.e., the 



some dynamical times 
is independent 



particular lattice: 

(i)a„iso(k,t>to)) 



(ei(k)-k)^ 



(40) 



where ei (k) is the eigenvector corresponding to the opti- 
cal branch, i.e., the one with maximal associated eigen- 
value. We plot this quantity in Fig. 21 Once again, we 
see that the bcc and fee lattices are very similar, while 
the breaking of isotropy of the sc lattice is much larger. 



V. CONCLUSIONS AND PERSPECTIVES 

In this technical paper we have explained step-by-step 
how to apply the Particle Linear Theory to a cubic Bra- 
vais lattice (sc, bcc and fee lattices) in a cubic simulation 
box. We use FFT techniques to speed-up the numeri- 
cal computations, which permits to compute the evolu- 
tion of the position of a large number of particles in a 
small computation time even with modest computer re- 
sources. We have illustrated the method computing the 
discreteness effects — in the linear regime — resulting 
from the evolution of continuous density field discretized 
using a perturbed bcc, fee and sc lattice. Attending to 
the tests we have performed, the bcc and fee discretiza- 
tions present less discreteness effects — in this regime 
— than the sc one, presenting small differences between 
them. They might be therefore better choices to set-up 
the IC in cosmological N-body simulations. 

As pointed-out in the introduction, an important mo- 
tivation of this work — and the reason for which we have 
actually developed it — is the study of the discreteness ef- 
fects in the highly-non linear (non-perturbative) regime. 
A way to estimate the discreteness effects in this regime 
consists in running a set of simulations set-up with dif- 
ferent Bravais cubic lattices [ij]. They lead to results 
which differ between them in a few per cent in the PS. 
From the IC and the final measured PS, a lot ingredients 
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enter in the game: the parameters of numerical integra- 
tion — strategy and accuracy in the computation of the 
force, smoothing, time-step — , finite-size effects, noise 
of the estimator, statistical fluctuations. . .It is impor- 
tant to have an analytic tool to check that the differences 
observed in the simulations corresponds actually to dis- 
creteness effects. The PLT plays this role in the linear 
regime of gravitational clustering (i.e. "small" k in the 
PS), which strongly suggests that the effects observed are 
actually discreteness ones in the whole range of k. 

We have described the method for an EdS universe. It 
is possible to genereralize the treatment for a flat back- 
ground model with cosmological constant — which is 
the currently most favored one — , without much extra- 
numerical cost. An implementation of the PLT for this 
kind of background model will be presented in a forth- 
coming paper. 

We have considered only cubic Bravais lattices. The 
method can be also applied to any Bravais lattice with 



the caveat that, in some cases, the simulations box could 
not always be a cube, i.e., the vectors R might not be 
translated into a cube using the transformation ([29|) . 
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APPENDIX A: EWALD SUM OF THE 
DYNAMICAL MATRIX 

The Ewald sum for the dynamical matrix is: 

V{R)=V^''\R)+V^''\R) (A1) 

with 



J 



P«(R^O) = -Gm^ 

n 

(5„,y 



GrnY^ 



(R-n-A)^(R-n-A)^ 



|R-n-A|3 



|R-n-A|2 

(R-n-A)^,(R-n-A)^ 
|R-n-A|5 



4q,3 

—^ exp(-a^ |R - n- Ap) 



erfc(a|R— n-A|) 
I 



2a 



(A2) 



exp(-a2|R - n-A^IR - n-A| 



and 



AnGr 



-D^ (R) = ^^9F E n^ ^-P f-S) ^- (k . R) fc.fc.. 



^- .^oH^I^ 



Note that the sum in Eq. (jA3[) is over all Fourier space 
and not only in the FBZ. The R = term is 

2?(R = 0) = - J2 2?(R-fn-A). (A4) 

In order to sum over a minimal number of vectors in real 
and Fourier space we take a « 2.067. 



In this particular case the functions [/„ (k, t) and Vn (k, t) 
can be calculated analytically and are: 



(A3) [/„(k,i)=5(k) 



^ \ "„ (k) ^ / ^ 

a.:,W I - j +a„(k) ( - 



4{^) 



T4(k,i)-5(k)io 



where 



to 



to 



^.Q-(k) .^x-a+(k)' 



(B3a) 
(B3b) 



APPENDIX B: SOLUTION OF THE MODE 
EQUATIONS 

We choose the solutions C/„(k, i) and y„(k, t) of the 
mode equation ([?I|) . without any loss of generality, sat- 
isfying 



a(k) 



C/„(k,to) = l, C/„(k,to) = 0, 
K(k,to)=0, K(k,io) = l. 
For an EdS universe the scale factor is [l|: 

2/3 



°'" ^ [t. 



(Bla) 
(Bib) 



(B2) 



and 



a„(k) 



a„(k) +Q„(k) 



v/1 + 24e„(k) - 1 



a+(k)--[v/l + 24e„(k) + l 



(B4) 

(B5a) 
(B5b) 



If £„(k) > the solution presents a power-law am- 
plification mode and a power-law decaying mode. If 
— 1/24 < e„(k) < 0, there are two decaying modes. Fi- 
nally, if £„(k) < —1/24, the solution is oscillatory and 



can be written as 



[/„(k,t) 



to 



7n(k)ln 



1 / t 



67„(k) V^o 



sm 



to 



7n(k)ln 



(B6a) 



to 



F„(k,t)=^!^(-| sm 



7n(k) Vto 



7n(k)ln(- 



where 



7„(k) = -v/|24e„(k) + 1|. 





(B6b) 



(B7) 



The evohition of the displacement field from any initial 
state u(R, to) is then given by the transformation 



u(R, t)^^Yl [^(^' ^)"(^' *o) + 2(k, t)ii(k, io) 



JkR 



(B8) 



where the matrix elements of the "evolution operators" 
V and Q are 



7'^,(k,t) =^C/„(k,t)(e„(k))^(e„(k)),, (B9a) 

n=l 
3 

Q^.(k,t) =^y„(k,t)(e„(k))^(e„(k)),. (B9b) 



n=l 



APPENDIX C: A BRIEF SUMMARY OF THE 
FFT TECHNIQUE 

Let us consider for sake of simplicity the one- 
dimensional FT / of the function /: 



3=0 



(CI) 



Because of the symmetries of the FT it is possible to 
divide the sum (|Cip (with N terms) into two sums with 
N/2 terms (this is called the Danielson-Lanczos lemma) : 



N/2-1 



A = E 



J27vjk/iN/2) 



2i 



j=0 

N/2-1 
^^2.k/N Y^ ^2.,k/(N/2)^^^^^^ (C2) 

J=0 



i.e., an "even" and "odd" term. Therefore, at this stage, 
it is possible to compute the even and odd sums at the 
same time, and then sum the result to obtain the desired 
FT. It involves a total of iV x {N/2) + 1 operations, in- 
stead of N"^ . For a number of particles which is a power of 
two, we can perform recursively the division (|C2p ln2 N 
times. Therefore the computation of the N terms fk in- 
volves only N \u2 N operations. This is called the Cooley- 
Tukey FFT algorithm. It exist other algorithms (which 
we will not describe here), which can use any number N 
of particles (and not only a power of two). 
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